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Abstract 

We present an ah initio approach for the computation of the magnetic sus- 
ceptibility x °f insulators. The approach is applied to compute x m diamond 
and in solid neon using density functional theory in the local density approxi- 
mation, obtaining good agreement with experimental data. In solid neon, we 

predict an observable dependence of x upon pressure. 
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The response of an extended system to a uniform external magnetic field is a fundamental 
property. This response can be used as a sensitive probe to the structural and electronic 
properties of materials, such as in the case of nuclear magnetic resonance spectroscopy. 
However, to our knowledge, the orbital magnetic susceptibility x °f rea l solids has not 
been computed from first principles. In this work we discuss an ab initio approach for 
the evaluation of x i n insulators within density functional theory (DFT). We applied our 
formalism to diamond and solid neon using the local density approximation (LDA) for the 
exchange and correlation energy. The agreement of our results with experimental data 
indicates that DFT-LDA describes correctly the magnetic response of these systems. 

The susceptibility x nas been evaluated in cubic semiconductors using empirical methods 
0. Exact expressions for x of a periodic solid in terms of Bloch eigenstates and eigenvalues, 
have been derived already in the sixties 0-§]]- However these approaches are rather involved 
and have not been applied to real materials. A more compact expression for x w &s recently 
given in Ref. ||, where it is applied to a model 2-dimensional system. Our approach for the 
computation of x in rea l systems is related to the one of Ref. |J. 

The paper is organized as follows. First we present the formalism for a generic single 
particle Hamiltonian. Then we justify the use DFT in the LDA in the computation of x-, 
and we discuss the accuracy and the limits of the additional use of the pseudopotential 
approximation. Finally, we apply our formalism to diamond and solid neon, studying the 
behavior of x as a function of the lattice constant. 

The magnetic susceptibility is defined as the second derivative of the total energy per 
unit volume E with respect to the macroscopic magnetic field B, i.e.: 

Xl3 ~ dBdBi [) 
where i and j are the Cartesian indexes. To simplify the notation in the following discussion, 
we consider a cubic system for which Xij — $ijX- 

Perturbation theory can be used to compute x- This is straightforward for a finite system. 
However, in an extended solid, the expectation values of the perturbative Hamiltonian on 
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delocalized eigenstates are not well-defined quantities for an uniform field. To avoid this 
problem we consider the response of the system to a magnetic field with a finite wavelength 
q = (g,0,0), i.e. B(ar) = b(0, 0, ^2 cos(gx)) = V x A with A(x) = b(0,V2sm(qx)/q,0). 
Defining 

in the limit of q — > 0, we obtain the macroscopic susceptibility x 0- 

Let us first consider a system described by a single particle Hamiltonian. If the coupling 
between B and the spin of the electron can be neglected, the perturbation to the Hamiltonian 
can be written as AH = + with 



(3) 



-Pyb 

c c q 

H m = _L A 2 = 1 s™ 2 (qx) b 2^ 

2c 2 c 2 q 2 



where atomic unit are used, p is the momentum operator, and c is the speed of light. 
For a periodic insulator we have: 
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where r/)^ and are the Bloch eigenstates and eigenvalues of the unperturbed Hamiltonian, 
Q is the volume of the unit cell, O and E are the sets of occupied and empty bands, and a 
factor of 2 for spin degeneracy is included. By inserting Eq. (0) in Eq. (f|), we get: 

2 r d 3 k 

xfa) = y ^t^ k + q ' k ) + ^ k - q ' k )] 

(5) 
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where N is the number of electrons per unit cell, 



and |iik,i) is the periodic part of the Bloch eigenstate (normalized in the unit cell). For 
q — > 0, the two terms on the right-hand-side (rhs) of Eq. (|5|) individually diverge, but x(q) 
remains finite. To show this, we use the f-sum rule: 

By inserting Eq. (0) in Eq. (f|) we obtain: 

2 /•d 3 k g(k + q,k)-20(k,k) + 3(k-q,k) 

= y ^ f • (8) 

Then 

X = lim X (g) = -|/^4" (k ' k ' )lk " k - (9) 
Similar conclusions have been obtained in Ref. ||. 

In our numerical evaluation of the macroscopic x we use Eq. (§) with a small but finite q. 
Note that Eq. (|5|) is not suitable to this approach. Indeed, in a practical application, both 
the integral in k space and the sum over all empty bands are replaced by finite sums. Under 
these conditions the f-sum rule, Eq. (0), is no longer exactly satisfied. Then for q — > the 
rhs of Eq. (^) will diverge as Af s c~ 2 q~ 2 , where Af s is the the error in the f-sum rule. This 
numerical instability does not occur in Eq. @ where every term is treated consistently. 

We computed x using DFT-LDA, i.e. we neglected any explicit dependence of the 
exchange-correlation functional on the current density. Ref. proposes an approximate 
functional for the exchange correlation energy E xc which depends also on the current. The 
current term in E xc influences the magnetic response in systems with a small electronic 
density. It is negligible in our case, since it yields a correction to x smaller than 2% at the 
electronic densities typical of the systems we are studying |||7J]. We also do not consider 
magnetic local field effects, which are negligible in non-magnetic materials ||. Finally, we 
note that the DFT Hamiltonian depends in a self-consistent way upon the electronic charge 
density. Thus, in general, to compute the second order variation in the total energy with 
respect to an external perturbation, one should take into account the linear variation of the 
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Hamiltonian induced by the linear variation of the charge Sp (see e.g. Ref. |IT|). However, 
if the perturbation is a magnetic field, 5p is zero by time reversal symmetry. Thus Eq. (||) 
is correct within DFT. 

In our present practical calculation we used the pseudopotential approach, in which 
only the valence electrons are considered. To discuss the validity of the pseudopotential 
approximation in the computation of x-, we divide the set of occupied bands O into the sets 
of core bands C and valence bands V. Then we have: 

X = Xc,e + Xv,£ = Xc~ Xc,v + Xv,e- (10) 

Here xc,£ is given by Eqs. (|6|) and (§) with the sum over the i and j indexes in Eq. (^) 
running over the the sets of core states, C, and of empty states, £, respectively. The other 
X with two indices are define in a similar way. xc is the magnetic susceptibility of the core 
electrons, which is not sensitive on the chemical environment and thus can be computed 
considering the isolated atoms, i.e.: 

Xc = Xc,e + Xc,v * ~ EE<*.VI*<>> ( n ) 
where we sum over the atoms in the unit cell, and ty\ are the core atomic wavefunctions of 
the atom /. Among the three terms in the rhs of Eq. (|TD|), Xv.£ is the only one accessible in a 
pseudopotential calculation; xc can be computed using an atomic code, but the evaluation 
°f Xcy requires the knowledge of both core and valence wavefunctions. Since xcy an d 
Xc are expected to be of the same order of magnitude, the pseudopotential approximation 
introduces an error of the order of xc by neglecting Xcy- This error is reasonably small 
only for elements in the first and second row of the Periodic Table, for which xc ^ X- F° r 
application of the present theory to heavier elements, all-electron calculations are needed. 
Finally, in our pseudopotential calculation, we replaced the operator — zV + k in Eq. (j^) 
with the velocity operator v£ = (d/dk)H^ where if£ is the pseudo-Hamiltonian flT3| . 

We computed x f° r isolated carbon (C) and neon (Ne) atoms, for solid Ne in the fee 
structure, and for solid C in the diamond structure. In the atomic phases we used the all- 
electron ground state wavefunctions to compute xc using Eq. (|TT|) . In Ne we also computed 



the atomic x using Eq. ( p|) with the sum over the index % running over all occupied states. 
In the solid phases, we evaluated xv,e using Eq. @ with a q = .037r/a, where a is the lattice 
constant of the cubic cell. The pseudopotentials were generated using the prescription of 



Ref. (T^|. In Ne we expanded the wavefunctions on a plane- wave basis set with a 120 Ry 
cutoff. We sampled the k space integrals with 10 special k-points in the irreducible Brillouin 
zone, and considered 400 empty states. In diamond we used a 60 Ry cutoff, 60 special k- 
points, and 300 empty states. We verified that with the above parameters the convergence 
error in the value of x is less than 0.2%. 

The results for Ne are shown in Table |. The atomic calculation is in good agreement 
with the experimental data. For the solid fee phase we report xv.e as a function of the lattice 
constant a. We note that xv,e reaches a plateau for a ~ a e Q) where Oq is the experimental 
equilibrium lattice constant. This indicates that for a > ag the interaction among Ne atoms 
is negligible. Moreover xv,e at a = ao is very close to the value of x computed for the 
isolated atom. This establishes, in the atomic limit, the correctness of our approach and the 
accuracy of the pseudopotential approximation. As a decreases, — xv,£ decreases. This can 
be understood by noting that for an isolated closed shell atom only a negative diamagnetic 
term contributes to x-, since the unperturbed Hamiltonian is spherically symmetric. As the 
Ne atoms get closer, spherical symmetry is broken and a positive paramagnetic term also 
contributes to x- F° r t ne sake of comparison with future experiments, we also report the 
theoretical pressure P as a function of a. Solid Ne at zero P is bonded by a weak van der 



Waals interaction, which is incorrectly biven by LDA [14]]. Thus for the larger a we do not 
expect to obtain accurate values for P. However we expect LDA to describe correctly the 
repulsive interaction between Ne atoms, which dominates P at smaller a. Note that at P=50 
GPa, — Xv,e is decreased by 16% with respect to its atomic value. 

The results for C are shown in Table Since C is not a closed shell atom, in the 
atomic case only xc is reported. In the diamond phase we report Xv,e a s a function of 
the lattice constant a. The computed pressure obtained from the LDA-DFT total energies 
is also shown. In the range of experimentally accessible pressures xv,s shows a negligible 
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dependence upon a. Both the values of xv,e at the experimental (<2q) and at the theoretical 
(oq) equilibrium lattice constant are in very good agreement with the experimental data. 

In conclusion we have presented a method to compute the magnetic response of real solids 
from first principles. We have shown that DFT-LDA reproduces the magnetic susceptibility 
X of diamond. In diamond x is found to be insensitive to the applied pressure whereas we 
predict an observable pressure dependence of x in solid Ne. 
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Institute for Basic Research in Science. Computer time was provided by the NSF at the 
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TABLES 

TABLE I. Magnetic susceptibility of atomic and solid fee Ne in units of 10 -6 cm 3 /mole. In the 
solid we considered different values of the lattice constant a. We indicate with a§ the experimental 
equilibrium lattice constant. The theoretical pressure P is also reported. 



-x 



-XC 



~XV,£ 



P (GPa) 



Atom (experiment) 
Atom (theory) 
Solid a =8.37au= ag 
Solid a =7.87au 
Solid a =7.37au 
Solid a =6.87au 
Solid a =6.37au 
Solid a =5.87au 
Solid a =5.37au 



7.2 
7.80 



.05 



7.79 
7.76 
7.64 
7.41 
7.14 
6.66 
6.04 



-2 
-2 
-1 
4 
15 
50 
151 
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TABLE II. Magnetic susceptibility of atomic C and of diamond in units of lCT 6 cm 3 /mole of 
C2. For the solid we considered different values of the lattice constant a. We indicate with a§ and 
Og the experimental and theoretical equilibrium lattice constants, respectively. The theoretical 
pressure P is also reported. 



-x 



~XC 



~XV,£ 



P (GPa) 



Solid (experiment) 
Atom (theory) 
Solid a =6.75au= <Zq 



Solid a =6.66au= Oq 



Solid a =6.55au 
Solid a =6.35au 
Solid a =6.15au 
Solid a =5.95au 



11.8 



0.32 



11.17 
11.23 
11.26 
11.23 
11.16 
11.09 



-17 


25 
85 
168 
283 
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